Macroscopic singlet states for gradient magnetometry 
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We present a method for measuring magnetic field gradients with macroscopic singlets states 
realized with ensembles of spin-j particles. While the singlet state is completely insensitive to 
homogeneous magnetic fields, the variance of its collective spin components is highly sensitive to 
field gradients. We compute the dynamics of this variance analytically for a chain of spins and also 
for an ensemble of particles with a given density distribution. We find an upper bound on how 
precisely the field gradient can be estimated from measured data. Our calculations can be used 
to carry out differential magnetometry with cold atomic ensembles using the multipartite singlet 
obtained via spin squeezing. 
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I. INTRODUCTION 

Realization of large coherent quantum systems is at 
the center of attention in quantum experiments with 
cold atoms [1, 2] and trapped ions [3]. Besides creat- 
ing large scale entanglement, it is also important to look 
for quantum information processing applications of the 
states created. Recently, a series of experiments have 
been carried out with cold atomic ensembles using spin 
squeezing [4] . This approach makes it possible to entan- 
gle 10 6 — 10 12 atoms with each other through making 
them interact with a light field and then measuring the 
light, realizing this way a quantum nondemolition (QND) 
measurement of one of the collective spin components 
[5, 6]. Spin squeezed states are useful for continuous vari- 
able quantum teleportation [7] and magnetometry [8-13]. 
In these experiments the atomic ensembles were almost 
completely polarized, which makes it possible to map the 
quantum state of these ensembles to the state of bosonic 
modes [14] and model with few variables even realistic 
dynamics including noise [15-18]. 

A basic scheme for magnetometry with an almost com- 
pletely polarized spin squeezed state works as follows. 
The total spin of the ensemble is rotated by a magnetic 
field perpendicular to it. The larger the field is, the 
larger the rotation is, which allows one to obtain the 
field strength by measuring a spin component perpendic- 
ular to the mean spin. So far, it looks as if the mean 
spin behaves like a clock arm and its position will tell us 
the value of the magnetic field exactly. However, at this 
point one has to remember that we have an ensemble of 
particles governed by quantum mechanics, and the uncer- 
tainty of the spin component perpendicular to the mean 
spin is never zero. Spin squeezing [19-22] can decrease 
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the uncertainty of one of the perpendicular components 
and this can be used to increase the precision of the mag- 
netometry [6]. 

Often the interesting quantity is not the absolute 
strength of the magnetic field, but its gradient and the 
effect of the homogenous field must be suppressed. For 
example, the magnetic field of Earth must be suppressed 
when the much smaller magnetic field around an elec- 
tric device or magnetic structure is measured [23]. The 
field gradient can be determined by differential magne- 
tometry, which can be carried out when two completely 
polarized atomic ensembles are used. In fact, the same 
light beam can pass through the two atomic ensembles, 
which can be used both for simultaneous spin squeezing 
of the two atomic ensembles and to carry out the differ- 
ential measurement [8, 9]. In general, singlet states of 
two large spins offer also the possibility for differential 
magnetometry [24]. 

It has recently been shown that interesting quantum 
states can be obtained even in unpolarized ensembles. In 
particular, if the uncertainties of the three collective an- 
gular momentum components are squeezed one after the 
other, then a multi-particle singlet state can be obtained 
[25, 26]. (For other approaches creating singlet states of 
cold atoms, see Ref. [27-29].) Singlets, as ground states 
of antiferromagnetic Heisenberg spin systems, have at- 
tracted considerable attention [30-34]. Such states are 
invariant under the action of homogenous magnetic fields. 
On the other hand, a magnetic field gradient rotates the 
spins at different locations differently, which leads to the 
gradual destruction of the singlet state. During this pro- 
cess, the variance of the collective spin components is 
increasing and this fact can be used to measure the field 
gradient. The advantage of this method is that only a sin- 
gle ensemble is used for differential magnetometry, rather 
than two ensembles, which leads to better resolution and 
also eases the experimental requirements on the method. 
Thus, ideas concerning the atomic magnetic microscope 
[10] and other techniques for sensing magnetic fields with 
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culations for a an atomic ensemble with a given density 
profile depicted in Fig. 1(b). In Sec. V, we consider the 
case of the singlet state mixed with noise and present re- 
sults for particles with a spin larger than i. Finally, we 
conclude the article in Sec. VI. 
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FIG. 1: (Color on-line) Multi-qubit singlet in a field gradient, 
(a) Equidistant chain of N particles, (b) Atomic ensemble 
with a Gaussian density profile. 



II. MULTI-PARTICLE SINGLET STATES 

In this section, we present an efficient description of 
multi-particle singlet states of N spin-^ particles. 

Pure multi-particle singlet states are eigenstates of J; 
with a eigenvalue for I = x,y, z. Here, the collective 
operators are defined as 
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where ji — t-&i with the Pauli spin matrices ai for / = 
x, y, z. Due to this, pure multi-particle singlet states are 
invariant under the unitary transformations 



Un{0) = CXp 



. ^n 



(2) 



where the angular momentum component along the n 
direction is 



atomic ensembles [11] can be applied also to gradient 
measurements. (Other methods for measuring the field 
gradient can be found in Refs. [13, 35].) The basic scheme 
for differential magnetometry with singlets is depicted in 
Fig. 1. 

Besides demonstrating the usefulness of multi-particle 
singlets for metrology, we also show that a certain type of 
quantum dynamics can easily be modeled for such states 
realized with an atomic ensemble. This is a surprise as 
quantum systems with millions of particles are typically 
difficult to model. 

While we mainly discuss spin-i particles, the spin 
squeezing procedure creating singlets and the differential 
metrology presented in our paper works also for spin-j 
particles for j > \- This is very important as looking 
for applications of ensembles of spin-j particles, without 
restricting the dynamics to the spin- i subspace is at the 
center of attention recently from the point of view of ex- 
periments and experimental proposals [25, 36-43], and 
also from the point of view of spin squeezing entangle- 
ment criteria [44-50]. 

Our paper is organized as follows. In Sec. II, we de- 
scribe the multi-particle singlet state of spin-^ particles 
and analyze its properties. In Sec. Ill, we calculate the 
dynamics of the variance of collective angular momen- 
tum components for such states under a magnetic field 
gradient for a singlet realized with a spin chain shown in 
Fig. 1(a). We also compute bounds on the precision of 
the gradient measurements. In Sec. IV, we present cal- 
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Such unitary transformations can be written as Ua(9) = 



,®N 



where u = exp (-| J2i nijiO) 



Mixed multi-particle singlet states are mixtures of pure 
multi-particle singlets. Hence, multi-particle singlets give 
zero for the expectation values of all moments of all col- 
lective angular momentum components 
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where I = x, y, z and m = 1, 2, ..., N. Mixed multipartite 
singlets are also invariant under the transformations of 
the type Eq. (2). In other words, they are in the zero 
subspace of the Hamiltonian 
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where k > is a constant. The dimension of this space 
is growing rapidly with N [53]. We have to identify the 
singlet created in the spin squeezing procedure in this 
space. 



A. Determining the singlet obtained in spin 
squeezing experiments 

In this subsection, we determine the multi-particle sin- 
glet created by spin squeezing procedures. Due to sym- 
metries of the setup, the state created is permutationally 
invariant. There are very many multi-particle singlets. 



Wc will now show that, on the other hand, there is a 
unique permutationally invariant singlet. 

Permutational invariance means that the quantum 
state g equals its permutationally invariant part 
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where Hk is a permutation operator and the summation 
is over all permutations. The singlet state realized by the 
squeezing procedure in an atomic ensemble is permuta- 
tionally invariant for the following reasons. First, It is 
created starting from the completely mixed state, which 
is permutationally invariant. Second, the measurement- 
feedback procedure to squeeze the collective variables in- 
volves only of collective, not individual variables. Hence, 
the dynamics is completely symmetric under exchange of 
particles [25]. 

Hence, we can state the following. 

Observation 1. For a given even number of particles 
N , there is a unique permutationally invariant singlet 
state. It can be expressed as 



(7) 
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t^o Tr ( e -^r 



Qs = (|*-> (*- 



l*->(*-l)PI, 
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where the operation (X)pi is defined in Eq. (6) and the 
two-particle singlet is 

l*-> = ^ (!+*>-£>.- !-*.+*>,). ( 9 ) 

where lis) are the eigenstates of j z . Proof. First 
we will show that there is a unique permutationally 
invariant singlet state. Such states have the follow- 
ing properties: (i) they are permutationally invariant 
and (ii) they are the eigenstates of Ji for I = x,y,z 
with eigenvalues 0. All permutationally invariant multi- 
particle states are uniquely characterized by the expecta- 
tion values (A®( N ~ n ) ® I 8 ™), where n = 0, 1, 2, ..., N- 1 
and A is a traceless single-particle operator [51]. More- 
over, as discussed before, all states for which J; = for 
I = x,y,z are invariant under the transformations of the 
type u® N , where u are unitary matrices acting on a sin- 
gle spin. Since any traceless A can be obtained from 
o~ z by unitaries, such a state can be uniquely character- 
ized by the N expectation values (a z ~ n ' <g> t® n ) , where 
n = 0, 1, 2, ..., N— 1. Knowing these expectation values is 
the same as knowing the expectation values of the pow- 
ers J", for n = 1,2,..., N. However, these expectation 
values are zero for all singlets as can be seen in Eq. (4). 
Thus, there is a single permutationally invariant singlet 
state, and Eq. (7) and Eq. (8) are indeed equal [52]. □ 



Let us interpret first the Eq. (7) formula. It de- 
notes a state that is a completely mixed state within the 
Yli(Jf) = subspace. It can also be written as 
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where \j, j z , a) denotes a state for which ^ Jf \j, j z , a) = 
j(j + 1) \J,3z,ot) , J z \j,j z ,a) = jz \j,jz,a) , a is used to 
label the degenerate eigenstates, and do is the degeneracy 
of the j = j z = eigenstate [53] . 

The alternative expression for the permutationally in- 
variant singlet is given in Eq. (8). It shows that the 
multi-particle singlet is an equal mixture of all tensor 
products of two-particle singlets [54] . One can even find 
that for an even N, the number of different permutations 
of such a singlet chain is [55] 

f(N) = (N- 1)!! = (N - 1)(7V - 3)(7V - 5) • ... (11) 

This fact has very important consequences for modeling 
quantum systems in such a state. While storing the den- 
sity matrix for a general quantum state of many particles 
is impossible, storing a representation of a product state 
or a state that is a product of few-particle units can be 
done efficiently. In the next sections we will explain how 
to compute quantum dynamics starting out from the per- 
mutationally invariant multi-particle singlet state. 



B. Calculating the reduced states of the singlet 
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In this subsection, we will calculate g[ cd , g 1 ^, and 
which are respectively the reduced one-particle, 
two-particle, and four-particle density matrices of the sin- 
glet state g s . Later, this will be needed when computing 
the time evolution of certain operators for the singlet 
state. 

We have to start from the decomposition Eq. (8). From 
that, we obtain the form 



P1234 



( a ^+^ 1*12) (*12|® 1*34) (*3 4 | 



+7 I ^12) (*i2 1 ® 7 + permutations 1,(12) 

where the second term has altogether 3, while the third 
term has 6 different permutations [including those ap- 
pearing in Eq. (12)]. The term multiplied by j3 is ob- 
tained when the particles 1 to 4 are in a product of two 
singlets states. The number of times that this special or- 
der is reached is given by the number of ways to distribute 
the remaining (N — 4) particles in (N — 4)/2 pairs, that 
is, (N — 4 — 1)!!. This has to be divided by the number 
of all distributions, that is, (N — 1)!!. Hence we obtain 



P- 



(N 



1)!! 



1 



(iV-1)!! (JV-l)(AT-3)' 



(13) 



The term multiplied by 7 is obtained when the parti- 
cles 1 and 2 are in a singlet state, but particles 3 and 
4 are not. As above, the number of times the remaining 
N — 2 particles can be distributed in pairs is (N — 2 — 1 ) ! ! . 
However, we have to subtract here the number of distri- 
butions where particles 3 and 4 are in a singlet state as 
well, that is, (N — 4 — 1)!! as shown above. Again, we 
have to divide this by the total number of distributions, 
(N — 1)!!, arriving at 
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(JV-1)!! 
1 



(JV-1)!! 

1 



(A-l) (N-l)(N-3) 



(14) 



Finally, for the coefficient of the completely mixed com- 
ponent we have 



1 - 3/3 - 67. 



(15) 



This occurs when all the four particles are in a singlet 
state with particles outside this set. 

The two-spin reduced density matrix can be obtained 
from Eq. (12) by tracing out particles 3 and 4 as 



Plf = Tr 34 (p 



d 
1234 
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where 



Ps = P + 7 = 



+ (l- ft )i (16) 



(17) 



N-V 

which can also be obtained from combinatorial calcu- 
lations similar to the ones carried out for /3 as (N — 
3)!!/(A— 1)!!. Finally, this leads to the trivial single-spin 
reduced state 

i 
2' 



P? d 



(18) 



III. GRADIENT MAGNETOMETRY WITH A 
MULTI-PARTICLE SPIN CHAIN 

In this section, we consider N spin-^ particles in a per- 
mutationally invariant singlet state of Eq. (8), where the 
particles are confined to a one-dimensional array. This 
could be prepared, for instance, with a Bose-Einstein 
condensate in an optical lattice driven to the so-called 
Mott-insulator state [61, 62]. 

We will calculate the effect of the magnetic field gra- 
dient on the singlet. In particular, we will calculate how 
the variance of the collective angular momentum compo- 
nents increases with the application of the field gradient. 
We will also calculate how precisely the field gradient can 
be estimated from the measured data. 

In our calculations, we consider a one-dimensional ar- 
ray along the ^-direction, such that the positions of the 
particles are given by 



where n = 1,2, .., N. 
tant chain where 



A particular example is the equidis- 



z n = (n — 1) ■ d + zo, 



(20) 



d is the distance between the particles and zq is an off- 
set. This situation is depicted in Fig. 1. We collect the 
positions on the z-axis in the vector z£- 

We can write the field at the atoms, situated along 
x — y = 0, as 



B{z 1 0,0) = B + zB 1 + O(z 2 ) 



(21) 



where we will neglect the terms of order two or higher. 
We will consider B = B • (0, 0, 1) and B 1 =B 1 - (0, 0, 1). 
For this configuration, due to the Maxwell equations, for 
the case of no currents or changing electric fields, we have 



divB = 0, 
curlB = 0. 



(22) 



This implies Yli= x y z^^i/dl = and dBi/dm — 
dB m /dl = for I ^ m. Thus, the spatial derivatives 
of the field components are not independent from each 
other. However, in the case of a linear chain only the 
derivative along the chain has an influence on the quan- 
tum dynamics of the atoms. Similar statement holds for 
a quasi one-dimensonal atomic ensemble, which is typi- 
cally the case if we consider an elongated trap. 

The Hamiltonian corresponding to the effect of a ho- 
mogenous magnetic field in the z-direction is 



N 



(23) 



where 7 is the gyromagnetic ratio. It gives rise to a time 
evolution given in Eq. (2) with ft = z, where z is the 
unit vector pointing in the z-direction. As we have dis- 
cussed before, multi-particle singlets are invariant under 
the transformations of the type Eq. (2). On the other 
hand, the singlet is not invariant under the quantum dy- 
namics generated by a magnetic field gradient B\ de- 
scribed by the Hamiltonian 



He 



N 



(») 



(24) 



where z c n is the position on the z-axis of the spin n and 
B\ is the field gradient along the z direction. Introducing 
a characteristic length L, the Hamiltonian Eq. (24) can 
be rewritten as 



H c 



WL 



iV 

E 

71=1 



(f>« 



(25) 



where Wj, = ^B\L. 

For instance, one may choose L = d in the case of the 
equidistant chain. Introducing the normalized Hamilto- 
nian H' r = — and 



{x n ,y n ,z n ) = (0,0,*£), 



(19) 



e = wit, 



(26) 



the time evolution operator becomes 



U G (Q) 



H' 

exp[-i— p-0]. 



(27) 



This formalism expresses the fact that our setup mea- 
sures the field gradient times the time. 

In order to use the singlet for differential magnetome- 
try, we need to find an observable that changes with 0. 
We investigate powers of collective operators J™ since 
these are relatively easy to measure experimentally. It is 
clear that I = z is not a good choice because [J z , H' G ] = 0, 
so this observable does not change with 0. Therefore, we 
start by investigating whether the first or the second mo- 
ment of J x might be a suitable candidate. 



A. Calculating (J x ) and (jf) 

We compute the dynamics of (J x ) and (J 2 ) starting 
from a multipartite singlet taking advantage of the fact 
that it is the mixture of all the permutations of tensor 
products of two-particle singlets as can be seen in Eq. (8). 
We will work in the Hciscnberg picture, thus all opera- 
tors will be given as a function of 0, while expectation 
values will be computed for the initial state g s . Hence 
the time evolution of the operator J x is given with the 
time-dependent single-spin operators as 



^(©) = E^ ) ( Q )- 



(28) 



In the Heisenberg picture, the time dependence of an 
operator A is given as 

A(@) = cxp (»§©) -4cxp (-&Q) , (29) 

Hence, for the time dependence of the single-particle op- 
erator ji 1 we obtain 



Jx V^V — c nJ x b njy 



e ' 



(30) 



where we introduced the notation 



c n = cos (y 1 ©) and s n = sin (y ©) ■ (31) 

The quantitites c„ and s n are the cosine and sine of the 
phases picked up by the n th particle. The expectation 
value of J x is therefore given by 



( J x )l N (©) 






=nOl n) ) 



^<4 n V ) = ° ( 32 ) 



because the single-particle reduced state of the singlet g s 
is the completely mixed state, cf. Eq. (18). Note that 
the reduced state p\ ed is the same for any n since g s is 
pcrmutationally invariant. In analogy, it can be shown 
that ( J y )s N (0) = 0. Therefore, a measurement of Ji is 



not suitable for estimating for I — x, y, z. We continue 
by investigating whether a measurement of J 2 is useful. 
We can write J 2 in the Heisenberg picture [cf. Eq. (28)] 

as a sum over two variables. Knowing that (Xq ) 2 = — , 
we can write the expectation value of the second moment 
of J x as a sum of a constant and two-body correlations 

MS 



(4 2 )s™(0) 



N , 



£ (4 ni) 4 n2) ) ft 

"1^12 
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In the following, we will partly drop the index z£ and 

sometimes also from (J X ) S N (0) when there is no risk 
of confusion. Since the singlet state is permutationally 
invariant, the correlation term on the right-hand side of 
Eq. (33) can be rewritten with the two-particle reduced 
state of the singlet as 

£ <4 rai) 4" 2) > es 

= E ((^ - «4 L, )(<w. C0 - * m 4 2) )W- 

(34) 

In Sec. II B, we obtained g\ e ^, the reduced two-spin state 
of the singlet [cf. Eq. (16)]. Direct calculation shows that 



(jk®3i)o& = 



4(iV-l) 



Ski, 



(35) 



where k,l = x,y, z and Ski is 1 if the two indices are equal, 
otherwise it is 0. Substituting Eq. (35) into Eq. (34), we 
obtain the sum of the two-body correlations as 



£ (4"°4 n2) ) es = 



ni/n 2 



4(JV- 1) 



where 






SnSmJ-t 



h, (36) 



(37) 



and therefore 



(4 2 )s Z "(©) 



Nh 2 



N{N - l)' 



(38) 



Since this is a non-trivial function of 0, a measurement 
of J x can be used to estimate the magnetic field gradient. 
Note that Eq. (38) contains a sum over two variables 
such that the two variables are not allowed to be equal. 
For practical purposes it is more useful to rewrite this 
with independent sums that require less computational 
effort as 

h = ^2(c n c m + s n s m ) - ^(4 + 4) = (C 2 + S 2 - N), 

n,m n 

(39) 



where we defined the sums 
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(40) 



and used that c 2 , -f 
cedure can be used 
the fourth moment 
below [56]. Insertin 
the following. 

Observation 2. 
tation value of the 
a singlet of a chain 
z^ , is given by 



s„ = 1. A similar subtraction pro- 
in a more complicated calculation for 
of an angular momentum component 
g Eq. (39) into Eq. (38), we can state 

The time dependence of the expec- 
second moment of J x , starting from 
of N particles with the z-coordinates 



(J')s N (e) 



Nh 2 
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N(N- 1) 



[N-C 2 -S 2 }\, 



(41) 



where C and S are defined in Eq. (40). 

Due to the symmetries of the setup, the dynamics of 
the variance of the y component of the angular momen- 
tum is the same as the dynamics of the variance of the x 
component, 



(^ 2 )s(e) = <j 2 ) s (0), 



(42) 



while (J 2 ) S (Q) = as mentioned before. 

From Sec. II we know that the singlet at t = is in- 
variant under the influence of any homogenous magnetic 
fields. For t > this is not true any more. However, 
the singlet evolving under the influence of Hq given in 
Eq. (25) still remains invariant under the transformation 
U- z {0) = axp(-*40) since [H G , J z ] = 0. 



1. Equidistant chain 

For the particular example of the equidistant chain, it 
is convenient to rewrite (J 2 ) s using the identity 



cos 



L 



-0 



(43) 



which from the Eqs. (38) and (39), and with the z- 
coordinates z^ e given in Eq. (20) leads to 



(4 2 >s(6) = 



Nh 2 



{i-^ny[E-([— ]*)-*]} 



(44) 
with the choice L = d. 

The variance of J x starts from zero as follows from the 
properties of the singlet state. It can also be seen from 
Eq. (44) because each of the N 2 terms in the sum is equal 
to 1 at = 0. Then the variance grows up to around the 




FIG. 2: (Color on-line) The dynamics of (J|)s/{Jj)wn (solid) 
and (AO)^ 1 (dashed) as the function of & for an equidistant 
chain of TV = 8 particles. (Jj)wn is defined in Eq. (45). 



level to that of the completely mixed state (white noise) 
1 
W 



(Jl) 



N, 



To be more precise, for even N 

Nh 2 



(J 2 ) s (& = n) 



1 



JV — 1 



(45) 



(46) 



which is very close to (J 2 ) wn for large N. The reason for 
this is that ^ n m cos ( [n — m\Kj =0 because when TV is 
even, the number of terms where n — to is even (such that 
cos (\n — m]ir) = +1) is equal to the number of terms 
where n — m is odd (such that cos (In — m]7r) = —1). 
At 6 = 2-7T, ( J 2 ) s returns to because again cos (\n — 
to]27t) = 1 for all n and to since n — to is an integer. 
Hence for the equidistant chain, (J 2 ) S (Q) is a periodic 
function with a period time 



T 



2tt 
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(47) 



If there is not a complete revival then the particles are not 
arranged in an array such that the intcrparticlc distance 
is uniform over the chain. This phenomenon can be used 
to characterize chains of atoms from the point of view of 
the uniformity of the distribution of the atoms. 

In Fig. 2, we plot the dynamics of (Jj) s (6) for N = 8 
spin-i particles. Note that the increase of (J 2 ) s is not 
mononotic for < 7r, it is rather oscillating. This is due 
to the fact that the atoms are arranged in a lattice. In 
the case of a continuous density distribution, there is not 
such an oscillation, as will be shown in Sec. IV. 



B. Calculating the precision of estimating 6 

From the calculation of (J 2 ) S (Q) it is clear that a mea- 



surement of J 2 gives information about 9, hence it can 
be used to estimate an unknown value of 0. In this 
subsection, we will calculate the precision A6 of the es- 
timation, for the singlet state realized with a chain of 
spins with the z-coordinates z£, undergoing a quantum 
dynamics due to a magnetic field gradient. 

The finite precision in reconstructing 9 comes from 
the fact that J 2 can only be measured with some uncer- 
tainty. The precision of the reconstruction of 9 based on 
measuring J 2 can be obtained as 



(A©) 2 



(AJ 2 



2\2 



\de(J* 



(48) 



where (AA) 2 = (A 2 ) — (A) 2 for an observable A. In order 
to calculate the expression Eq. (48), we need to know 
(J 2 ) S (Q), which we have just obtained, and (J 4 ) s (9), 
which wc arc going to calculate now. 

Let us consider again the one dimensional chain of par- 
ticles. For this, the Hamiltonian is given in Eq. (25). As 
in Sec. Ill A, we will work again in the Hciscnbcrg pic- 
ture. The time evolution of the operator J 4 is 



■jfi 2 , we arrive at 



(4) 
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+h 2 



N + 3N{N - 1)] 
3(N - 2 " 



E <*e 

^(ni,n 2 ,n3,rc.4) 



E (4 ni) 4 n2) k 

ni^n2 
(ni) Y (n 2 ) v {n 3 ) y(ni)\ 

A e A e A e /a> 



(52) 



The expectation value of the sum of the two-body corre- 
lations is given in Eq. (36). Next, we calculate the expec- 
tation values for the four-body correlations. For that, we 
use the reduced four-particle density matrix presented in 
Eq. (12). For the reduced four-particle matrix we obtain 



(-.(I) ,-(2) „-(3) ,-(4) 

\Jx J x J x J x I g 1 ^ 



aw)< 



16(JV-l)(AT-3)' 
ft 4 3 

16(AT-l)(AT-3)' 



(53) 



^(@) = (E^ n) (0) 



(49) 



The expectation value of this can be rewritten as a sum 
over four variables as 



and 



/o(l) 7 -(2),-(3) 7 -(4)\ 
\Jx Jx Jy Jy I c 



ft 4 



16(JV-l)(iV-3) 



(54) 



<•£>. 



ni,n2,n3,rn 



The singlet state is invariant under [/® w for any lo- 
cal unitary U. This leads to (j x j y j y j y ) e rcd = 
(X^X^X^X^)^. (50) yWjWjWjW)^ =0[60]. Since' the singlet S 3 per- 



Again, we leave out the indices zfi and 9 for simplicity 
at the moment. We can rewrite Eq. (50) with sums in 
which the variables of the summation are not allowed to 
be equal as 

(™i)\4\ 



(j 4 ) s = £<(4 ni) ) 



-3 J2 ((4 ni T(^" 2) ) 



»1^12 



+4 y: ((4 ni) ) s (4 Ba) )> 



»1^12 



+6 Yl ((4 ni) )ww 3) )> 



^(n 1 ,n 2 ,n 3 ) 

E 

^{n 1 ,n 2 ,n 3 ,n 4 ,) 



I V^ /J[(ni) -gfa) <g(na) gin*) 



e /& 



(51) 



where Yl-t/j 7 m denotes summation over the indices j, k, I 
such that none of them is equal to another one. Based 
on Eq. (51) and by making use of the fact that (Xq ) 2 — 



H II 

^1234 

mutationally invariant, with this we know all the expec- 
tation values with four-body correlations with j x and j y . 
Taking everything into account, we arrive at the expres- 
sion 



(J 4 ) s = |J3iV 2 



2N - 
3 



67V -8 
N -I J 

h 



(7V-l)(7V-3) 



(55) 



where we defined 

h= E [ckClC n Cm 

^(k,l,m,n) 



SkSlS n S m + 2ckClS n S m ]. (56) 



The term J4 involves a sum over four variables, which 
implies a large computational effort for large systems. 
Similarly to how we have done with I 2 in Eq. (39), we 
can write I4 as 

h = E \- c k c i c nC m + s k sis n s m + 2c k cis n s m -P}. (57) 

k,l,n,m 

In P we will write all the terms that appear in ^ fe ( 
but do not appear in J2^(k imn)- Let us study the first 



term of the right-hand side of Eq. (57). We can rewrite 
it as 

/ J [ckQCnCm + SkSlSnS m + 2CfeC(S„S m ] 
k,l,n,m 

= C 4 + S 4 + 2C 2 S 2 , (58) 

where C and S are defined in Eq. (40). Next, we have to 
determine P. After determining it, we have to eliminate 
the sums containing conditions that multiple indices are 
unequal in a similar manner by replacing the sum with 
another one without such a condition and substracting 
the difference. Finally, we arrive at an expression with 
single-index sums only. Using these results, we obtain 
the following [56]. 

Observation 3. The time dependence of the expec- 
tation value of the fourth moment of J x , starting from a 
singlet of a chain of N particles with the z-coordinates 
z^ , is given by 



(4) z s N (Q) 
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+ 



2N- 
3 



6_/y__8 

X-l J 
h 



(59) 



(X-l)(iV-3) 



where Ii and I4 are defined in the Eqs. (37) and (56) and 
can be rewritten as 



2 = 


= X^ + X^-N, 


4 = 


= -X"i,o + -^0,1 + 2 ^i,o^o,i 




—6X4 — 6X ,4 — 12X 2 ,2 + 3X 2 + 3X 02 




+8X3,0X1,0 + 8X0,3X0,1 + 4X :1 + 8X2,1X0,1 




+8X1,2X^0 + 2X 2 ,oX ,2 — 6X 2 ,qX 1 — 6X ,2X j 



-2X2, oX 1 — 2Xq, 2X X q — 8X1,1X1,0X0,1, 



where 



N 



Xk,i = J2 c « s « 



(60) 



(61) 



ra=l 



with c n and s n as defined in Eq. (31). 

Alternatively, it is possible to write the term I4 in a 
more compact form as [56] 



X <{ 2(X - 3) - 4X(X - 2) /i(a) 

4 



-X J 



/i(«) 



X 



h(2a) 



-2X 2 Re[/ 1 2 (a)A(2a)*]}, (62) 



where 



A (a) 



V^ e iaz l 

N ^ 
k 



and a 



e 

T 



(63) 



which is also easier to compute. 

Observation 2 and 3 make it possible to calculate the 
precision of the estimation A6 based on Eq. (48). It is 
possible to calculate analytically the precision for = 0. 
For that, we determined the zeroth and first order terms 
of the Taylor expansion of Eq. (48) using the expansion 
of sine and cosine up to second order. Hence, we can 
state the following. 

Observation 4. The maximal precision of estimat- 
ing AG for an equidistant chain of N spin-^ particles is 
characterized by 



(A6)- 2 (6 = 0) = 



X 2 + X 3 
12 



4 



N 



X- 1 



(64) 



where a = L^/(N 2 — 1)/12 is the standard deviation of 
the equidistant chain, and L = d. 

A more general result on the precision at = will 
be presented in Observation 7 below. 

In Fig. 2, we plot the dynamics of (AO)^ 1 for X = 8 
for the equidistant chain. As can be seen, the precision 
is maximal at O — 0. Then, it decreases, however, this 
decrease is not monotonic and the precision is oscillating. 
In particular, based on Eq. (48) one can see that the pre- 
cision is zero when the tangent of (Jj) s (0) is horizontal, 

i.e., a e (4 2 ) s (e) = o. 



IV. CONTINUOUS DENSITY PROFILE 

In this section, we work out the formulas describing the 
case of a one-dimensional continuous density profile. We 
present the dynamics of the second and fourth moments 
of the collective angular momentum components for this 
case. 

In the case of a spin chain, the particles were placed in 
the fixed positions z^. Now, while we still consider the 
case when the particles are localized in certain positions, 
the distribution of X particles is given by a distribution 
function /jy, where 

/jv(zi, -22, — , z N )dz 1 dz 2 -dz N = f N (z N )dz N (65) 

is the probability that particle 1 is between z\ and z\ + 
dzi, particle 2 is between zi and Z2 + dz2, etc. Without 
loss of generality, //v can be considered as invariant under 
the permutation of any two particles [59]. As before, 
we compute the average (J 2 ) s in order to estimate the 
magnetic field gradient and also (J^) s in order to estimate 
the uncertainty A0. For a general distribution function 
/jv(2jv), they are given by 



(^>s /N (0)= dz N f N (z N )(J^(Q) 



(66) 



with (Jj)f N (0) from Eq. (41), and 

(4) f s N (©) = / dz N f N (z N )(4f s - (0), (67) 



with (J%)l N (O) from Eq. (59). Here and in the following, 
we substitute z£ with Zn in the expressions for fixed 
particle positions. 

In order to compute particular examples, we need some 
properties of the reduced M-particle distribution func- 
tions fm (M < N). These can be obtained from /jv 
as 



Im{zm) = / dz M +i-dz N f N (z N ), 



(68) 



where we used the shorthand notation zm = 
(zi, z 2l ...,zm) t as above. Note that due to the invari- 
ance under permutations of /jy, it does not matter which 
of the N — M particles are integrated over. Further, Jm 
is permutationally invariant as well. 



A. Distributions with Dirac 8 functions 

For example, the distribution function of the chain con- 
sidered in the previous Section can be written as [56] 



1 N 



jV! ^ J-l " V j k i 

' ^(k 1 ,k 2 ....,k N )j = l 



Its reduced distribution functions are given by 

M 



(69) 



fM N (zM) = 



(N- m)\ 

AH 



e nn^-^ 



^(k 1 ,k 2 ,...,k M )j=l 



(70) 



The pcrmutational invariance of /m simplifies many 
calculations. For instance, the average of I 2 from Eq. (37) 
for a distribution /at(zat) becomes 



(h 



dz N I 2 (z N ) 



(71) 



= J^ / dz N f N (z N )[c n c m + s n s m ] 

/ / QZ n <AZ m j2\Z n7 Z m )[C n C m T" S n S m \ 



n^m l 



- N{N - 1) / dz 2 f 2 (z 2 )[c lC2 + Sl s 2 }. 

In analogy, we obtain the average of I4 from Eq. (56) for 

Jn(zn), 



(h) f * 



dz N f N {z N )h{z N ) 



(72) 



df 4 / 4 (z 4 ) 



(iV-4)! 

X [C1C2C3C4 + S1S2S3S4 + 2C1C2S3S4]. 



In the following, we will make the assumption that the 
system is a gas of particles that are uncorrelated in space, 
which means that the distribution function can be writ- 
ten as a product of single-particle distribution functions, 



F N {z N ) = K=ih{z n ) 



(73) 



Further, we will assume that f\ is a smooth function 
which does not contain Dirac delta functions, i.e., points 
with infinitely high density. 



B. Independently and smoothly distributed 
particles 

Let us compute the expectation values of J 2 and J 4 as 
we did in the previous Section, for a product distribution 
function of the form given in Eq. (73). Considering the 
Eqs. (38) and (55), it becomes clear that we only need to 
compute {I 2 )f« and (J 4 )^. Since / 2 p (z 2 ) = /i(z 1 )/ 1 (z 2 ) 
in the uncorrelated case, we obtain from Eq. (71) that 



(I 2 ) f « = N(N~1)(C 2 + S 2 ), 



where 



and 



C 



s 



dzi/i(zi)cos ( y6 



dzi/i(zi)sin 







(74) 



(75) 



(76) 



With these, we have all the ingredients to determine the 
dynamics of (J 2 ) s for the gas from Eq. (38). 

Observation 5. For an ensemble of particles with 
a product distribution function f^ from Eq. (73), we 
obtain the following dynamics for the expectation value 
of the second moment of J x : 



(Jl)l»{&) 



Nh 2 



1 - c 2 - s 2 



(77) 



where C and S are defined in Eq. (75) and Eq. (76). Here 
the averaging is over the density distribution fi(z\). 
In analogy, we obtain from Eq. (72) that 



(h) f " = 



N\ 



(JV-4)! 



(C 2 + S 2 ) 2 



(78) 



Inserting (I 2 )* N and (I^)* 1 * into Eq. (55), we obtain the 
dynamics of (Jx)s- 

Observation 6. For an ensemble of particles with 
a product distribution function f^ from Eq. (73), we 
obtain the following dynamics for the expectation value 
for the fourth moment of J x 



(4 



4\/K 



Nh 4 
To" 



3N - 2 - (6N - 8){C 2 + S 2 ) 
+3{N-2)(C 2 + S 2 ) 2 \, (79) 
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where C and S are defined in the Eqs. (75,76). 

Note that due to the product structure of /JJ-, this 
expression is much simpler than the expression for the 
chain given in Eq. (59). These formulas make it possible 
already to calculate the precision of the phase estimation. 



1. Gaussian density profile 

One of the most common density profiles is the Gaus- 
sian profile. As an example, we will now calculate the 
dynamics for this case explicitly. The Gaussian density 
profile is given as 



/r uss (^i) 



(*i-y) 2 
zC a? 2 



V27TCT 2 



(80) 



where zq is the coordinate of the point with the highest 
density and a is the width of the profile. Substituting it 
into Eq. (77), we obtain 



C 2 + S 2 = e—& 



(81) 



and hence 



<^(e) = ^ 



1 — TT©' 

1 — e i2 



Note that the maximum of ( Jj) s (6) is exactly the value 
for white noise, cf. Eq. (45). Thus, in this case there is no 
overshooting, the noise does not become larger than that 
of the white noise, as was the case with the equidistant 
chain in Eq. (41). 

In Fig. 3, we plot the dynamics of ( J 2 ) s for N = 10 5 
spin- j particles. The density profile is Gaussian with 
a = L. 



330 



220 



© 
< 



110 




0.024 



0.016 



0.008 °V 



0.05 0.1 





FIG. 3: (Color on-line) The dynamics of (Jj) s / '(Ji)wn (solid) 
and (AB)^ 1 (dashed) as the function of for N — 10 5 par- 
ticles with a product distribution function of Gaussian pro- 
files [Eq. (80)] for a = L and (Ae)^ 1 calculated. We also 
present (AS)^ 1 calculated based on the Gaussianity assump- 
tion Eq. (108) discussed in Sec. VB (dotted). (J^) wn is de- 
fined in Eq. (45). 



(82) For the derivation, see the Supplementary Material [56]. 
For any product distribution function f^, (AO)~ 2 (0 = 



0) = Njjz since cov(zi,z 2 ) = in this case. Hence 
the maximal sensitivity is a simple function of N and 
a. Correlations can increase or decrease the sensitiv- 
ity compared to the uncorrclated case. For instance, 
the correlations present in the equidistant chain increase 
(A9)- 2 (9 = 0) by a factor ^-, cf. Eq. (64). 



FURTHER CONSIDERATIONS 



2. Analytic formula for the maximum precision for any 

particle distribution, including correlated particle 

distributions 

We will now give an analytic formula for the maximum 
precision at O = for the case of any particle distribu- 
tion. 

Observation 7. The maximal precision of estimating 
A6 is characterized by 



(AG); 2 (G = Q) = N 






cov(zi,z 2 ) 



L 2 



(83) 



where 



a 2 = 



cov(zi,z 2 ) 



dzfi(z)(z-(z)) 2 , 

dzfi(z)z, 

dzubfiizi, zn){zi - ( 2 l))( 2 2 - {Z2})- 



(84) 



A. The influence of noise 

So far, we considered only perfect singlet states. In 
practice, the multipartite singlet state cannot be realized 
perfectly. In this section, we will consider the case of 
starting from imperfect singlet states. 

A realistic method to model the imperfections is by 
introducing local decoherence channels for each qubit as 



) (g) = {l-q)g + q 



1 

2' 



(85) 



where < q < 1. Each spin is mixed with a certain 
amount of white noise locally. The above single-qubit de- 
coherence is given with the Kraus operators as €q (g) = 
(1 -q)g + qY,i K iQ K l where K = {t,a x ,o- y ,a z }/2. 

We assume that all these decoherence channels act in 
parallel for all spins. Thus, we obtain the multi-spin 
decoherence 



£ q wn) (Q) 



4 wn ' 1) (^) o 4 wn ' 2) (e) 



°4 Wn,iV) (0)- (86) 
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In Eq. (86), o indicates function composition (fog)(x) = 
f(g(x)) and the number in the superscript indicates to 
which qubit the function is applied. We consider an 
atomic ensemble with a continuous distribution. Next 
we outline briefly the derivation of the noisy dynamics 
for this case. For simplicity, we assume that the decoher- 
ence acts after the unitary dynamics. 

The moments (J 2 ) and (J x ) given in Eqs. (33) and 
(52), respectively, are affected through the two-body and 
four-body correlations via the formulas 



E 

ni/ri2 



and 



<4 B1 



X 



tllo 



(l-?) 2 



E <*t i: 



X, 



(«2) 



711^712 



(87) 



\ A e A e A e A e )e. 



I Y"( rai ) v( n 2) v-(™3) y-( n4 )\ 

\ A e A e A e A e lo* 



^(ni, 712,713,714) 

= a-^) 4 E 

/(il,i2,n3,«4) 



Hence, for the atomic cloud we obtain (cf. Observation 
5 and 6) 



(J 



2\/X 



Nh 2 



x/s,g 



l-(l-q) z (C z + S z ) (89) 



and 



(Ji)& --■- ^\3N-2-(l-q) 2 (6N-8)(C 2 + S 2 ) 



'x/s,g 
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+ (l-q) i 3(N-2)(C 2 + S 2 ) 2 \. (90) 



Note that for q = 1 Eqs. (89) and (90) reduce to the 
values corresponding to a global white noise. Eq. (89) 
can be rewritten as 



(4 2 kwn(<z) = (l-qf(J 2 x ) s + ll-(l-q) 2 }(J 2 x ) 



i/wn- 



(91) 



Let us consider a Gaussian density profile described in 
Sec. IV B 1 for which we have C 2 + S 2 given in Eq. (81). 
Substituting Eq. (81) into Eqs. (89) and (90) we can com- 
pute the influence of the noise for the second and fourth 
moments of the angular momentum components. The 
precision of the reconstruction of O based on measuring 
J 2 can be obtained as 



vanish at = 0. In contrast, the numerator is a positive 
number. Hence, even for very small amount of noise we 
have 



(A0)- 2 n (0 = 0) = 0. 



(93) 



Note that this is due to the fact that we chose a noise 
state which is invariant under the O dependent transfor- 
mation Ug from Eq. (27). The precision remains close to 
zero until the added noise is much larger than the noise 
of the singlet. 

(ii) The other limit is the case of the large O. Let us 
define wn such that for > wn the singlet evolved 
into a state that, based on the second and fourth mo- 
ments of the angular momentum coordinates, is like the 
completely mixed state. That is, for for > wn we 
have 



(Jl) 






(Jl) 

(Jt) 



(94) 



The reason for this is that for large enough the quan- 
tity C 2 + S 2 as given in Eq. (81) is close to zero and in 
Eqs. (89) and (90) only the constant terms corresponding 
to the moments of the white noise remain. Hence, based 
on Eq. (92) for > wn , we have 



(A0)-L(0)«(l- (? ) 4 (A0)- 2 (0). 



(95) 



In Fig. 4, we calculated the precision (A0) s (0) 
as a function of the noise for q = 0.01,0.05 and 
0.1 for a Gaussian ensemble. For these values, 
(A0)" 2 (0)/(A0)-2 n (0) = 1.04,1.23, and 1.52, respec- 
tively. 

The noise we considered in this section is a particular 
type of noise. Other local noise realistic in some scenarios 
is 



4 x) (Q) = (l-q)g + q\ + 1 1 ) x ( + 



(96) 



Modeling the effect of such a noise is more complicated 
than modeling the local white noise even if the noise is 
acting after the unitary dynamics in our model. However, 
for a realistic modeling of the noise Eqs. (96) we have to 
assume that the noise acts in parallel with the unitary 
dynamics because the state of the particles affected by 
the noise given in Eq. (96) is influenced by the magnetic 
field and has a dependence. In contrast, the particles 
affected by white noise were insensitive to the magnetic 
field. 



(A0) 



2 

s.wn 



Ph2 



(4 4 kwn(0) - (4 2 )s,w n(0) 

(l- q r\d e (J 2 ) s 



(92) 



Let us see now the behavior of the precision Eq. (92) 
under noise in the limiting cases, for q > 0. 

(i) At = 0, the denominator of the right-hand side of 
Eq. (92) is zero. This can seen from Eq. (44), from which 
it follows that d(j(J 2 ) s is a sum of sine expressions, which 



B. Spin-j particles 

So far we have discussed the case of singlet states of 
j = h particles. In this section we will study the singlet 
of spin- j particles. We will find that the dynamics of (J 2 ) 
for N spin-j particles is the same as the dynamics of (J 2 ) 
for N spin- i particles, when in both cases we normalize 
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FIG. 4: (Color on-line) The dynamics of (Jf )s/{J% )w» (solid) 
and the precision (AO)^ 1 for the noiseless case (dashed) and 
for a noise with q = 0.01,0.05 and 0.1 (from top to down, 
dotted) as the function of 6 for N = 10 5 particles. (J|) wn 
is defined in Eq. (45). The density profile is Gaussian and 
a = L. 



Note that {j[ l) j\ 2) ) s 



Ui Ji )s,j f° r n ^ n' due to 



the permutational invariance of the state. Hence we ar- 
rive at 



(jl n) il n \j 



(Ji 



(100) 



N(N - 1)' 
and we can rewrite Eq. (35) for particles with j > i as 



(jk ® jl) r>*,3 



(Jl) 



N(N - 1) 



Vd , 



(101) 



red,.J 



where g™ 2 ' J is now the reduced state of the spin- j singlet 
[60]. 

The time-evolved single particle operators j x (&) = 
X(Q) have the same form as in the j = | case, given 
in Eq. (30). Therefore, we can perform the calculation 
for any j in analogy to the derivation for the spin- 2 case 
of Sec. III. A. Starting from 



(■4 2 (0)> 



s,j — \J x )wn,j 



ni/n 2 



(ni) v (n 2 )\ 



'X 



<-) 



/sji 



(102) 



the variance with that of the white noise. Moreover, we 
find that using a certain Gaussian assumption in order 
to estimate ( J x ) s , we obtain the same dynamics even for 
the precision. 

When considering the spin-j case, one could think 
about using the ideas of Observation 1 to obtain the 
multi-particle singlet of spin-j particles as an equal mix- 
ture of all the permutations of tensor products of two- 
particle singlets. However, the state obtained this way 
does not equal Eq. (7). This is due to the fact that for 
j > 2 there are several permutationally invariant SU(2) 
singlets [63]. Hence, for j > \ another method is needed. 

For the dynamics of the variance of the collective angu- 
lar momentum components, we need again the variance 
for the completely mixed state (white noise), which is 
obtained for N spin-j particles as 



(Ji 



m±R# 



I 2 ) 



(Ji 



where (J%) 



(97) 



is defined in 



Note that (J 2 
Eq. (45). 

Let us denote by ji the angular momentum compo- 
nents of a single spin-j particle. Since for each spin the 
dynamics is invariant under a coordinate transformation, 



,-(™)\2 



( n )\2 



from (jrT + (jf T + (iz? = j(j + l)^ 2 we obtain 



<£w l) ) 2 x ) 



(•"i/wiji 



(98) 



for I = x,y,z and for all 0. Moreover, from the require- 
ment that ( Jf) s ,j = 0, we obtain 



(Ji 



(£tfi B) ) a )»J + N(N - lXjfjfk; = 0. (99) 



we arrive at 



where 



[<-£>(e)kj = ",•[<•£> (Q)L, 



(Jl) 



(J 2 ,) 



ri(i + i)- 



(103) 



(104) 



Therefore, Observation 2 generalizes to spin-j particles 
as follows. 

Observation 8. The dynamics of the variance of J x 
for a chain of spin-j particles with the z- coordinates z^ 

is 



(J!)Z (©) = «j 



h 2 N 



l 



(105) 



where C and S are defined in Eq. (40). 

For an ensemble of particles with a density profile \(z), 
we can in analogy generalize Observation 5 to any j as 
follows. 

Observation 9. For an ensemble of spin-j particles 
with a product distribution function f^ from Eq. (73), we 
obtain the following dynamics for the expectation value of 
the second moment of J x : 



(j 2 Jhe) = Kj ^ 



i-cr-S' 



(106) 



where C and S are defined in Eq. (75) and Eq. (76). Here 
the averaging is over the density distribution fi(z\). 
For the Gaussian density profile Eq. (80), we arrive at 



(J 2 x )(Q) = ^h 2 (l-e-^ 2 



(107) 
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later fits very well the exact dynamics. 

Invoking the Gaussian assumption (J*) = 3(J 2 ) 2 for 
the precision calculations for N spin-j particles we obtain 
the same result as for N spin-i particles 



(A9)„-(e) = (Ae),i(e). 



(109) 



With these considerations, even with the calculations 
with the local noise q carried out for j = ■= fit the calcu- 
lations for j > i carried out for the same level of noise. 



VI. DISCUSSION AND CONCLUSION 



FIG. 5: (Color on-line) The dynamics of (Jf ) S /{J% )«.» (solid) 
and (AO)^" 1 (dashed) as the function of for N — 6 particles 
in an equidistant chain. We also present (A0)~ calculated 
based on the Gaussianity assumption Eq. (108) discussed in 
Sec. VB (dotted). 



The calculation of (J x )(<3) for the singlet of spin-j par- 
ticles seems to be much more complicated than for spin-^ 
particles. It is possible to avoid calculating the fourth- 
order moment by using the assumption that when J x is 
measured, the probability of the measurement outcomes 
follow a Gaussian curve, with the zero outcome being the 
most probable. For such Gaussian probability distribu- 
tions, the higher order moments can be obtained from 
second order ones [58]. In particular, for our case 



(J 



HJl) 



(108) 



Equation (108) leads to (AJj) 2 w 2(jJ) 2 , keeping in 
mind that (J x ) = 0. Such a Gaussianity assumption is 
expected to work for later times, as at 9 = the vari- 
ance of J x is zero, and in the beginning only few of the 
eigenstates of J x are populated. Later, however, many 
eigenstates of J x are populated and a continuous approx- 
imation of the discrete spectrum is appropriate. Note 
that the Gaussianity of the probability distribution is a 
notion completely independent from the Gaussianity of 
the density profile of the cold gas. 

One can substitute Eq. (108) into the formula for 
(AO) 2 given in Eq. (48) which can be used for any j. 
This clearly simplifies the calculations. In Fig. 3, we also 
plot the dynamics based on the approximation using the 
Gaussian assumption on the correlations Eq. (108) for 
the j — i-case as an illustration. The bound obtained 
this way diverges at t = 0, however, it is very close to 
the true value for > 0.025. We also carried out a cal- 
culation for j > |. In Fig. 5, you can see the comparison 
for the case of a chain of six atoms. The precision based 
on the Gaussian assumption diverges at t = 0, however, 



We presented calculations for a many-particle singlet in 
a magnetic field gradient for an ensemble of distinguish- 
able, well localized spin-i particles. We have shown that 
multi-particle singlets created in cold atomic ensembles 
can be used for differential magnetometry. The magnetic 
field gradient can be estimated by measuring the vari- 
ance of one of the components of the collective angular 
momentum. We calculated the dependence of this vari- 
ance on the field gradient and the measurement time. We 
also calculated the precision of the estimation of the field 
gradient. We have also considered admixtures of white 
noise and discussed the extension of the results for gen- 
eral spin-j particles. Our work opens up the possibilities 
for experiments with unplorarized ensembles. 

In the future, it would be interesting to find bounds 
for the precision of the field gradient measurements us- 
ing the theory of the quantum Fisher information [64]. 
It would also be important to consider the case of parti- 
cles that are not well localized, in which case the spatial 
degree of freedom cannot be easily separated from the 
internal degrees of freedom for the quantum dynamics 
we considered [59]. For room temperature experiments 
with atomic ensembles, the atoms must be well localized, 
however, singlets can also be realized with Bosc-Einstein 
condensates, for which all particles are delocalized. For 
such systems, the noise is very different from the case of 
distinguishable particles. 
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Magnetic microscope for differential magnetometry (Supplementary Material) 

I. ADDITIONAL CALCULATIONS FOR OBTAINING (J*> 8 (6) FROM OBSERVATION 3 

The general expression of (J^) s (0) for N particles sitting at the positions z%, collected in a vector z^, has been 
obtained in Eq. (59). Here and in the following we will label with zfi a vector with N elements z\, ...,z c N . In order to 
being able to compute this for large N, we need to simplify 



-*4 ■ -*cccc ~r -*ssss T" ^^ccssi 



(SI) 



where 



-*cccc ■ — / CkC-lC-mCn i 

^(k,l,m,n) 

-*ssss • / SkSlSrnSn-) 

-'cess * = / C]^CiS m S ni 

^(k,l,m,n) 



(S2) 



where c^ = cos(-y-0) and Sfc = sin(^0). We will show two ways of doing this. Firstly, as stated in Eq. (60), one can 



rewrite it as 



-^1,0 + ^o,i + 2 ^i,o^o,i 

— (6^2,0-^1,0 + 6^0,2^0,1 + 2J 2 ,o^o,i + 2X ,2^i,o + 8^1,1^1,0-^0,1) 

+8-^3,0-^1,0 + 3X 20 + 8X 0i3 X 0yl + 3X 02 + 8X^1X0,1 + 2X2,0^0,2 + 8X lj2 ^i,o + &X\ 4 

-(6X4,o + 6Xo,4 + 12X 2 , 2 ), 



(S3) 



where 



N 



(S4) 



The proof is presented in Section I A below. Secondly, as stated in Eq. (62), one may also rewrite it more compactly 
as 



I 4 =n\2(N-3)-4N(N-2) /i(o) + N 3 fr(a) +N /i(2a) - 2JV 2 Re[/ 1 2 (a)/i(2a)* 



(S5) 



where 



A(«) = ^E e<a * J and a = r- 



(S6) 



The proof is presented in Section IB below. 



A. Proof of Eq. (S3) 
Let us concentrate on the last term in Eq. (SI). We can write using a shorthand notation 

^ccss I / 

\k,l,m,n 

- E 

- E 

^(k,l,m— n) 

E E . E E 

y^(k—m,l,n) ^(k—n,rn,n) ^(k,l—m,n) ^(k,l—n,m) 

E 

7^(fc— Z,?ti— n) 

E - ,E 

E - E 

^(^k—l—m,n) ^(k—l—n,m) 

E - E 



^ c k cis m s n . (S7) 



fc=z= 



m—ny 



Here ^ (k = m,l = n) means that the summation is such that k = m,l = n and k ^ I. Eq. (S7) can be rewritten 
after simple considerations as 



-fccss = 2J C kClS m S n - ^ ( c l s l s m + CkClS^ + 4c k S k ClS m ) 
k,l,m,n ^(k,l,m) 

- Y ( c k s i + ScfcSfcQS; + 2c 2 k s k ci + 2c k cisf) - Y cj^. (S8) 

Each term in Eq. (S8) corresponds to a line in Eq. (S7). Then, we can rewrite the terms in Eq. (S8) still containing 
the conditions "not equal" with terms without such conditions as follows 

E c l siSm = [12- 12 - 12 - 12 - 12 ) c k siSm 

=£(k,l,m) \k,l,m ^L{k=l,m) =£{k,l=m) =£(k=m,l) k=m=nl 

= 2_^i C k S l S m 



k.l. 



2(E <%*ksi- Y fei ) - (E c ^ - E c ' s fc) 






E2 2 
c k s ky 



Y c fc s ' s ™ - E( 2c ' sfcs * + °k s ^ + 2 E c ' s '' ( S9 ) 



k.l.m kj k 



where we used that we have 



/] akh = y] akbi - y^ aubk (S10) 

^(k,l) k,l k 

for any real numbers a k and b k . Analogously, one finds that 

Y^ c kCis 2 m = ^ c kCis 2 m -J2(2c k cisf + c 2 k sf) + 2^2 (%sl, (Sll) 

^{k,l,m) k,l,m k,l k 

^ C k S k ClS m = ^ c kSkClS m - }^(cls k Sl + C k S k ClSi + C k Cisf) + 2 ^ c l s l- ( Sl2 ) 

^(k,l,m) k,l,m k,l k 

Substituting Eqs. (S9), (Sll), and (S12) into Eq. (S8), and using again Eq. (S10) for the remaining terms of two 
non-equal indices, we arrive at 



and 



,m,n 

^ c his m + J2{24 s kSi + c 2 k s 2 ) - 2 Y2 c l s 



k.i. 



kA,m k.l 



^2 c k cis 2 n + J2(2c k cis 2 + c 2 k s 2 ) -2j2c 2 k s 

k,l.ni k,l k 



-4 ^2 c k s k cis m + 4^(cfeS fc sj + c k s k cisi + c k cisf) -8^0^^ 

k,l.m k,l k 

- ^2{c 2 k s 2 + 2c k s k c 4 si + 2c 2 k s k s l + 2c k cis 2 ) + 7 ^ c 2 k s\ 

k,l k 

-£44. ( S13 ) 

In Eq. (S13), the first four lines correspond to the first line in Eq. (S8), and the remaining two lines to the second 
line. This can be simplified by combining terms that appear more than once as 

-'cess — / J CktySmSn 
k,l,m,n 

- 22 ( C k s l s m+C k CiS 2 n +Ac k S k ClS m ) 
k,l,m 

+ J2 (ctsksi(2 + 4 - 2) + cls 2 (1 + 1 - 1) + CfcQS?(2 + 4 - 2) + c fc s fe qs,(4 - 2)) 

k,l 

+ J2c 2 k si(-2-2-8 + 7-l), (S14) 

k 

which finally yields 

^ccss = ^2 C kClS m S n ~ ^2 { c t s l s m + C^iS^ + Ac k S k ClS m ) 
k,l,m,n fc,Z,m 

+ ^2 {^ c l s kSi + c 2 k s 2 + Ac k cisf + 2c k s k cisi) -6^2 c k s k- ( S15 ) 



k.l 



The formula for / cccc can be obtained from the formula for / ccss [Eq. (S15)], by replacing s k by c k and combining 
terms that appear more than once as 



/cccc= J2 C kClC m C n - J2 C l C l°™ (1 + 1 + 4) + J2 C k°l ( 4 + 4 ) + ^ c fcC?(l + 2)- &J2 C k 
k,l,m,n fc,/,m k,l k,l k 

= ^2 c k cic m Cn - Q^2c 2 k cic m + £(8cfq + 3c 2 k cf) - 6^4- (S16) 



k.l.m.n k.l 



Similarly, the formula for 7 SSSS is obtained as 

/ ssss = 53 s k sis m s n - 6 y^ s\sis m + y^(8a|gj + 34Q -6^4- (S17) 

k,l,m,n k,l,m k,l k 

Combining the results of the Eqs. (S15), (S16), and (S17), we obtain 



/ CfcQCm^n ■+" SkSiS m S n ~r ZCkCiSmSn 

k,l,m,n 

~ X! ( 6c fc C * C ™ + &S 2 k SlS m + 2c 2 k SlS m + 2c fc QS, 2 „ + 8CfeSfcCjS m ) 

+ ^ (8cgci + 3c^c^ + 8*g*i + 344* + 8<^a fc *j + 244 + 8c k cisf + ±c k s k asi) 

k,l 

-6^(4 + 4 + 244). (sis) 

This is equivalent to Eq. (S3). rj 

B. Proof of Eq. (S5) 

Using the continuous distribution formalism one can write J4 from Eq. (SI) as 

Ji = (M - A)\ / d ^ 4 -fi N (^ 4 ) ( C 1 C 2 C 3C4 + S1S2S3S4 + 2ciC2S 3 S 4 J , (S19) 

where j\ N (£4) is the reduced 4-body correlation function for the chain, cf. Eq. (72) of the main text. It is computed 

as f A N (£4) = J dz§ • • • dzwf N N (zn) from the permutationally invariant N-vaiiaic probability density of N particles 
with the ^-coordinates z^, which is given by 

N N 

& <*») = TV- d n s 0* - **°w) = w e n s (* - 4) ■ ( s2 °) 

7res N fc=i 7 t(fc 1 ,fc 2 ,...,fc N )i=i 

Here, SV is the permutation group of N particles and the first sum runs over all permutations ir from that group. In 
the second sum the kj indices ranges from 1 to N with the restriction that the indices be different, and z k are the 
locations of the particles on the chain. 

For any permutationally invariant probability density /jv one can show that 

/ Az A fAc l c 2 czCA L + S1S2S3S4 + 2C1C2S3S4J = / di4/ 4 cos 6 = / dz^f^e 1 = e (S21) 

holds. The second equality holds because the sine expression occuring in the exponent is an odd function under the 
exchange of Z\ + zj, and z 2 + £4. In this way, we can express I4, with the help of characteristic functions. In particular, 
the multivariate characteristic of the multivariate probability density /at(£at) is 

f N (a N ) = ^££=1°**) = [dz N f N (z N )e^™=i akZ ». (S22) 

As can be easily checked, /jv(cE/v) has the following properties: (i) /jv is permutationally invariant if /jv is, (ii) 

/am [o ; i,a2, ..., Qtjv-ii 0] ) = /at-i(q?jv-i), (hi) /at(Ojv) = 1; where Ojv is a vector where all entries are equal to 0, and 

(iv) Jn(-^n) = /jv(cfiv)*- Property (ii) follows from the fact that J dz/v/)v(£v) = /jv-iO^v— 1) and property (hi) 
follows from the normalization of /n{zn)- 



Comparing the Equations (56), (S21), and (S22), we observe that 







(N - )l iN ( a ' _a ' a ' _a )' where a - L 



(S23) 



and where f% N has to be computed from the probability density of Eq. (S20). In general, the lower elements f^ (zm) 
(M < N) are given by 



nn^i) 



{N - M)\ 



M 

E IP (*-*%) 



(S24) 



^(k 1 ,k 2 ,...,k M )j = l 



Let us compute the characteristic function of / A ,™ , dropping from now on the upper index zfi in order to simplify the 
notation. We obtain the following recurrence relation 



/jVf(aM) = 



(TV -A/)! 

m 

(N-M)l 

m 

(N-M)\ 



E 



M „ „0 



'££i 



^(k 1 ,k 2 ,...,k M ) 



E « 

JV( e * aMZ1 >' 



£ 1 <*,•*£ 



M-l 

E E 

^(fei,fc 2 ,...,fe M _ 1 ) i=i 



^EllT^.^.^aM^ 



EM-l o \ 

j = i a i z i\ - 



Af-1 



N -M 



(N-M+1)W I (AT-M + 1)! f^ 

— — — - < JV/i (ajvf) /m-i(<3a/-i) - ^2 /m-i(«m-i + a^e;) > , 



'1^=1 a i z i ^,a M z l 



(S25) 



where e; is a vector of length M — 1 that has only one nonvanishing element (that is equal to 1) at the position I. It 
can be used to compute 1^ via Eq. (S23), leading to 



/ 4 (a, -a,a, -a) = |JV/i(a)/3(a,-a,a) -2/ 2 (a,-a) + / 3 (a, -2a, a)| 



(S26) 



We can apply the recurrence relation again for M = 3 in order to reduce the complexity of this expression. We obtain 

1 



/s(a3) 



{iV 2 /i(ai)/i(a 2 )/i(a3) + 2/ x (ai + a 2 + a 3 ) 



(AT- l)(iV-2) 

-^[/i( Q i)/i( a 2 + a 3 ) + /i(a 2 )/i(ai + a 3 ) + A(a3)/i(«i + a 2 )] | , 
which for the two cases of interest in Eq. (S26) yields 



/ 3 (a, -a, a) = 
f 3 (a, -2a, a) = 



1 



(TV- l)(JV-2) 
1 



N z 



/i(a) /i(a) - JV/i(2a)/r(a) - 2(7V - l)/^) 



7V 2 /i 2 (a)/r(2a)-27V /:(«) - TV A (2a) 



(N- 1)(N -2) 
Similarly, we obtain for M = 2 that 

/ 2 (ai,a 2 ) = |iV/i (ai) f x (a 2 ) - A («i + a 2 )j 

For the special case of interest a.\ = — a 2 occuring in Eq. (S26) this reduces to 

1 



+ 2 



f 2 (a, -a) 



JV-1 



JV 



/i(«) 



(S27) 



(S28) 



Finally 

/ 4 (a, -a, a, -a) 



(TV- 1) (AT -2) (AT -3) 



2(iV-3) -4AT(AT-2) /i(a) 



TV d 



/i(a) +iV A (2a) 



~2A^ 2 Re 



[/?(«)A(2a)*]}, 



which due to the identity (S23) is equivalent to Eq. (S5) for /i(a) = -^ J2k e laZh [computed with the Eqs. (S22) and 
(S24)] with a= §. n 

II. ADDITIONAL CALCULATIONS FOR OBTAINING (A6) s " 2 |e=o FROM OBSERVATION 7 

We will show that for 9 — > 0, the inverse variance of the estimation of 9 is given by 



N 
L 2 



(A9) s z |e=o = t^ [<J 2 - cov( Zl ,z 2 )] 



(S29) 



where 



a 2 = /dzi/i(2i)(zi-^i)) 2 , 

(^l) = / dzi/i(zi)zi, 
cov(zi,z 2 ) = / dzidz 2 /2(^i,2 2 )(zi - (z 1 ))(z 2 - (z 2 )) = {ziz 2 ) - (z 1 )(z 2 ). 



Proof. We estimate the uncertainty from the error propagation formula 

(AJ 2 ) S 



(A9) 



\de(J*) s 



(S30) 



(S31) 



cf. Eq. (48) of the main text, for general continuous density profiles. The quantities which occur are (Jj) s (9), 
<9e(^)s(9) and (J 4 ) s (9). In order to get the desired limit, we need to expand them around 9 = 0. Let us start with 
(J 2 ) s (9). For fixed particle positions zn, we obtain 



(4 2 >r(e) 



Nh 2 



l 



N(N- 1) 



n? 



n^rn 



- / n *"m 



L 



(-) 



from the Eqs. (37,38) and (43). Averaging this over a general pcrmutationally independent density profile /n(zn), 
we obtain 



(J 2 x )t N (Q) 



Nh 2 



1 — / dzidz; 



1/2(21, z 2 ) cos f 






9 



Nh 2 



(1-J2). 



Expanding the cosine in the integral we arrive at 

1 



la Ml 



2L 2 



• d Zl dz 2 f 2 ( Zl , z 2 )(zi - z 2 ) 2 9 2 + (9(9 4 ) 
1 " ^2 ( / d«i/i(«i)«i " / dz 1 dz 2 / 2 (z 1 , z 2 ) Zl z 2 ) 9 2 + 0(9 4 ) 



-l(a 2 -cov(z 1 ,z 2 ))9 2 + 0(9 4 ) 



(S32) 



where the last line is obtained by adding and subtracting a term (zi) 2 . We also used that due to the permutational 
invariance (z 2 ) = (z 2 ) and (zi) = (z 2 ) hold. This leads to 



NH 



(Jl)l N (©) « ^ (> - cov(z! , z 2 )) e 2 + 0(9 4 ) 



and 



5e(J 2 ) s /N (9) 



Aft 2 
2I 2 " 



a 2 -cov(z 1 ,z 2 ) 9 + 0(9 3 ) 



(S33) 



(S34) 



Let us now consider the expansion of the term (J 4 ) s (0). Again for fixed positions zn, we have [cf. Eq. (67)] 

(4)f N (6) 3iV 2 -2iV W - 4 ^ (Zn-Zn 



4 y^ (z n ~ z rn \ 6 1 
— r > cos " H -T -i 4 , 

1 ) f-f V L J 16(N-l)(N-3 ' 



ft 4 16 8(7V-1)£^ V L 7 16(JV-l)(JV-3)' 

with J 4 from Eq. (SI) above. Note that in contrast to this equation, the particles positions are labelled by zn instead 
of z£ because we have to average the expression over fpj. This leads to 

(J 4 )H©) _ 3N 2 ~2N JV(3JV-4) /",,,, v (z x -z 2 



h 4 16 

3N(N -2) 



/ d2;idz 2 / 2 (z 1 ,z 2 ) cos f ^— 2 6 
/ d 4 z/ 4 (z 4 )( C1C2C3C4 + S1S2S3S4 + 2C1C2S3S4 



16 

3JV 2 - 27V iV(3JV-4) - 3N(N-2) f 
= —16 8 /2 + 16 /4 ' (S35) 

where we used again the permutational invariance of /at. We need to expand the expression J 4 . Using the first 
equality from Eq. (S21) and expanding the occurring cosine as before one obtains that 

h « 1 - 2^2 / d 4 z/ 4 (zi, z 2 , z 3 , z 4 )(zi + z 2 - z 3 - z 4 ) 2 e 2 + o(e 4 ) 

= 1 - -^[ct 2 - cov(zi, z 2 )]6 2 + 0(6 4 ). (S36) 

Inserting the expansions of I2 from Eq. (S32) and of 7 4 from Eq. (S36) into Eq. (S35) leads to 

(■g)fr(8) , 3^-2iV_7V(3iV-4 )( 1 

ft 4 16 8 V L 2L v ' ;J 



+ m * 6 2) (l - ^k 2 - covfr, z 2 )]6 2 ) + 0(9 4 ) 



iV 
4L 



2 [a 2 -cov(z 1 ,z 2 )]e 2 + 0(6 4 ). (S37) 



Now we have all the necessary ingredients to prove the claim. Indeed, inserting the Eqs. (S33,S34,S37) into Eq. (S31) 
we obtain 

, A ^-2 „ ggV - cov^, z 2 )] 2 6 2 + Q(6 4 ) N , unrftJ . 

fjVi " - cov(zi,z 2 )J0 2 + O(0 4 ) ^ 

which proves Eq. (S29). rj 



